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Abstract 

By a new Monte Carlo algorithm we evaluate the sidedness prob- 
ability p n of a planar Poisson-Voronoi cell in the range 3 < n < 1600. 
The algorithm is developed on the basis of earlier theoretical work; 
it exploits, in particular, the known asymptotic behavior of p n as 
n — > oo. Our p n values all have between four and six significant dig- 
its. Accurate n dependent averages, second moments, and variances 
are obtained for the cell area and the cell perimeter. The numerical 
large n behavior of these quantities is analyzed in terms of asymptotic 
power series in n _1 . Snapshots are shown of typical occurrences of ex- 
tremely rare events implicating cells of up to n = 1600 sides embedded 
in an ordinary Poisson-Voronoi diagram. We reveal and discuss the 
characteristic features of such many-sided cells and their immediate 
environment. Their relevance for observable properties is stressed. 

Keywords: planar Voronoi cells, Monte Calo algorithm, large sided- 
ness 
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1 Introduction 



A Voronoi diagram partitions space into convex cells constructed around a 
set of point-like 'seeds' or 'particles', in such a way that each point of space 
is in the cell of the particle to which it is closest. When the particles are 
distributed randomly and uniformly, the partitioning is called a Poisson- 
Voronoi diagram or a Random Voronoi Froth. 

Voronoi cells play a role throughout science and engineering and are also 
of interest to mathematicians. Applications include cellular structures that 
either arise spontaneously in nature (e.g. in biological cellular structures, in 
soap froths, or in granular materials) or are employed as a tool of analysis 
(e.g. to identify lattice defects in simulations of melting crystals). Many 
references are given in Ref. pQ and in the encyclopedic monograph on tessel- 
lations by Okabe et al. 

The simplest Voronoi diagrams are of the Poisson type. It is important, 
therefore, that the properties of Poisson- Voronoi diagrams be understood as 
well as possible. Here we pursue, by means of a new Monte Carlo method, 
earlier investigations [fl [3j H] on such diagrams in the Euclidean plane M 2 . 

The most prominent statistical property of the planar Poisson- Voronoi 
cell is its 'sidedness'. We denote by p n the probability that a cell is n-sided, for 
arbitrary integer n > 3. Other properties of interest include the average area 
of an n-sided cell and the average length of its perimeter; the statistics of the 
angles at the vertices; and correlations between neighboring cells. All these 
properties may be expressed as multiple integrals on the particle positions 
[51 [2J, but only a few of them can be calculated explicitly. In particular, 
no simple closed form expression for p n is known. An exact relation derived 
from Euler's theorem ensures that the average sidedness n = J2^=3 n Vn is 
equal to n = 6. 

It is known numerically that p n peaks at n = 6 and falls of rapidly for 
large n. Hayen and Quine [6] have numerically evaluated the integral for p 3 
with high accuracy. For n = 4, 5, . . . the values of p n stem only from Monte 
Carlo work. The most accurate reported values of p n are due to Calka [7] 
for 4 < n < 7 and to Brakke [S] for 8 < n < 16. One has pi$ ~ 10~ 8 and 
the largest sidedness observed in simulations by conventional algorithms is 
around n — 16. Drouffe and Itzykson [HUTU], as part of an effort to construct 
field theories on random lattices, developed an improved algorithm by which 
they estimated p n for n up to 50. Their results, however, have error bars that 
for n ^ 30 become of the same order as the p n themselves. Hence simulating 
many-sided Voronoi cells has remained a challenge. 

The interest of investigating Voronoi cells for asymptotically large n was 
stressed by Le Caer and Delannay [TT]. Analytic knowledge of the large 
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n behavior of p n , apart from the insight that it provides, also constrains 
the laws that describe the finite n behavior as observed in experiments and 
simulations. An example of this interplay between the regimes of finite and 
of asymptotic n is the theoretical explanation given in Ref. [I] of the failure 
of Aboav's law [12] for Poisson-Voronoi diagrams. 

The analytic study of p n in the limit n — > oo was taken up in Refs. [5] 
and [lj. It was shown there, among many other things, that asymptotically 

p n ~ Cp®, n - oo, (1.1) 

with C = 0.344 347... P] and 

Pn ~ 4vr 2 (2n)! • 1 ] 

In the present work we exploit this asymptotic knowledge. Going beyond 
Eq. (11.11) we write an equality that is exact for all n rather than merely 
asymptotic, namely 

p n = C nP $\ (1.3) 

whence necessarily lim^oo C n = C. In this work we focus on C n and show 
that it can be expressed as an average 

C n = (6e- v ), (1.4) 

where V is a known expression in the angular variables that describe the 
n-sided cell, and 6 is an indicator (i.e., equal to or to 1) imposing a 
geometric constraint on the set of angles. We will determine the prefactor C n 
in (11.31) by Monte Carlo evaluation of the right hand side of Eq. (jl.4p for finite 
n = 3,4, . . .. The Monte Carlo algorithm is new for this problem. Whereas 
all previously used methods become rapidly inefficient with increasing n, the 
performance of the algorithm presented here is, very roughly, independent of 
n. This makes it possible, in particular, to explore the structure of Voronoi 
cells in the hitherto inaccessible large-n regime. 

The remaining sections of this paper are the following. In Sec. [2] the al- 
gorithm is described. In Sec. [3] results are presented and discussed for the 
sidedness probability p n as well as for the averages and second moments of 
the cell perimeter and cell area. The asymptotic large-n behavior of these 
quantities is analyzed numerically. In Sec. H] we present and discuss charac- 
teristic pictures of many-sided Voronoi cells in an environment of ordinary 
cells. In Sec. Owe summarize our results. 

The algorithm requires the explicit expressions for V and G in Eq. (jl.4p . 
Finding these is a matter of considerable technical complexity; it is based on 
results of Ref. jT] and is the subject of Appendices IA1 and IB1 
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2 Monte Carlo algorithm 



2.1 Context 

Monte Carlo methods for generating Voronoi cells of Poisson distributed 
particles are discussed in the monograph by Okabe et al. [2J. One class 
of methods simply determines p n as the relative frequency of occurrence of 
n-sided cells. But since p n decreases to zero faster than exponentially for 
n > 12, the statistical precision goes down accordingly With such methods 
it is hardly possible to accumulate sufficient statistics for even single-digit 
precision as soon as n « 16. 

Another class of methods generates cells for a value of n fixed in advance. 
The first to have done so seem to have been Drouffe and Itzykson [9j. The 
method employed by Calka |7J is also in this class. These methods face the 
problem of attrition: a Monte Carlo generated geometrical object, in order 
to represent a valid n-sided cell, must satisfy certain geometrical constraints. 
The probability that an attempted generation satisfy the constraints again 
decreases rapidly with growing n. 

The present algorithm, which also fixes n in advance, completely solves 
the problem of attrition: the geometric constraints are satisfied with a prob- 
ability that tends to unity when n — > oo. In order to arrange things this 
way, a certain amount of rather technical rewriting of the initial problem is 
necessary. We have confined this rewriting to the Appendices. If one accepts 
its results, the method is easy to apply. 

2.2 Angular variables 

An n-sided Voronoi cell around a particle in the origin, as shown in Fig.fU is 
specified completely by its n vertex vectors Si, S 2 , . . . , S n . It may be specified 
alternatively by its n mid-point vectors, i.e. the projections Ri, R 2 , . . . , R n 
of the origin onto the sides. The explicit expression 0, [EH 13 EE] for p n as a 
multiple integral on the R m is given in Appendix [A] It has not, however, 
been possible to evaluate this integral analytically. By choosing other sets 
of variables of integration one may recast the original integral in numerous 
different forms, none of which is simple. For our purpose it is essential to use 
the angular variables that we will define now. 

Let $ m and \l/ m be the polar angles of R m and S m , respectively. Other 
angles relevant for this problem are defined in Fig.[TJ The r\ m = \P m+ i — ty m 
are the angles between two consecutive vertex vectors and the £ m = $ m — 
$ TO _i those between two consecutive vertex vectors; n-periodicity in the index 
m is understood. For fixed sets £ = {£ m } and r\ = {q m } one may still jointly 
rotate the set of vertex vectors with respect to the set of mid-point vectors: 
this modifies only the relative angles (3 m and 7 m between the two sets. We 
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may select any one of these relative angles and call it 'the' angle of rotation, 
since it will determine all others; we will select (3\. When for a generic (3\ 
we draw the cell boundary by clockwise constructing its successive segments, 
then after a turn of 2tt it appears not close onto itself but to spiral. A 'no- 
spiral condition' must therefore determine the appropriate value of the angle 
of rotation fit for which the cell boundary closes and which we will denote 
by Pi = /?*(£, 77). This condition reads 

G(^;&) = (2.1) 

where 



c 



2irG TT cos ^ 



n^r- (2-2) 

771=1 

One may note that Eq. (12.21) involves the (3 m and 7 m that are themselves 
determined by the solution fii = (3* of (12.1 1) . For an arbitrary pair (£,77) 
there need not exist a solution to Eq. (12.11) . In Appendix [B] we show that it 
has a solution, which moreover is unique, if and only if 



m—1 m—1 



max 

Km<n 



- Vi) + Zm ~ min < 7T, (2.3) 

^— ' J Km<n L ' 



1=1 ~ ~ 1=1 



which is a criterion expressed entirely in terms of the supposedly given sets 
£ and 77. 

After these preliminaries we return to (jl.4p . The symbol G in that expres- 
sion denotes the indicator function of the domain in (£,77) space where (12.31) 



is satisfied. Finally, the 'interaction' V in (11.41) is given explicitly in terms 
of the angular variables in Appendix |A] through a sequence of definitions, 
Eqs. f lA.lOp and flA.5l) - flA.8l) . that we will not display here. 



2.3 Algorithm for determining p n 

The sidedness probability p n is given by Eqs. (I1.2l) - (ll.4p . We determine it 
numerically by evaluating (Oe _v ) as follows. We fix the sidedness n, after 
which the simulation proceeds as follows. 

(i) We draw n — 1 random numbers uniformly distributed on [0, 1] and 
order them. After multiplication by 27r this gives [15] < ^1 < ^ < • • • < 
^ n -i < 2n. We set ^ n = 2n and choose 

Vm = ^m+l-^m, 771=1, ...,71-1, 

Vn = *1- (2.4) 

We next draw 2n — 1 random numbers, order them, and discard those of 
odd rank so that only n — 1 are left. After multiplication by 2ir this gives 
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Figure 1: Heavy line: the perimeter of the Voronoi cell around a particle in the 
origin O. Dashed and dotted lines connect the origin to the midpoints R m and 
vertices S m , respectively. The particles of the neighboring cells are located at 
2Ri,...,2R n . Right angles have been marked. The figure defines the sets of 
angles £ m , rj m , /3 m , and -y m . 

< $1 < <l 2 < . . . < $ n -i < 27T. We set 4> = and choose 

£ m = $ m -$ m _i, m = l,...,n-l, 

£ n = 2vr -$„_i. (2.5) 

(ii) We check if the pair of sets (£,rj) thus obtained satisfies Eq. ( 12. 3D . If 
so, then we know that there exists a ?y) which may be determined from 
Eq. (12.11) . hence = 1 and we proceed with (iii). If not, then it is impossible 
to satisfy Eq. (12.11) . we have = 0, and the attempt to generate an n-sided 
cell fails. We increase the attempt counter by one unit and return to (i). 

(iii) We solve /?*(£, 77) from Eq. (12.1 ft by a numerical iteration procedure 
which also yields the derivative G'(£,r); /?*) needed in the next step. 

(iv) We calculate the weight exp(— V) according to Eqs. (lA.lOj) and (]A.5j) - 
( 1A.8I) of Appendix [A] and add it to the accumulated weight. We increase the 
attempt counter by one unit and return to (i). 

(v) In the end the total accumulated weight is divided by the total number 
of attempts, including those that failed. The result is an estimate for p n . 

We remark that the successive cells generated by this procedure are all 
statistically independent. 
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2.4 Algorithm for n dependent averages 

The simulation method described above allows us to study arbitrary cell 
properties F(Ri, . . . , R n )- Writing (F) n for the average of F subject to a 
given sidedness n we have 

(F) - ^ Qe ~ V > (2 6) 

{F)n ~ (0e-v> ' (2 - 6) 

Here the numerator, which generalizes (11.41) . has an insertion Ip that derives 
from F by a radial integration. We recall that the average (...), defined in 
(1A.15j) . applies to quantities that depend exclusively on the angular variables. 
To find Ip from F, we let R av denote the average of the R m . Upon setting 
R m = RsvPm we may express the p m entirely in terms of the angular variables 
(see Appendix |A]) . Then, if F is of dimension d F , it may be factorized into 
a radial and an angular part according to 

F(R 1 ,...,R n ) = (C/4A)^ /2 F(e,r ? ), (2.7) 

where we show explicitly the areal particle density A which had heretofore 
been scaled away [16]. When (12.71) is integrated over the radial scale R av , an 
extra factor appears as compared to the same operation for p n and we find 

T(n + hd„) j in ^ 
I F = V — 2 F> W- d r /2 F , 2.8 
Y(n) 

where V denotes the Gamma function and where we abbreviated 

W = 4A7r(l + n~ l V) (2.9) 

with V given by (IA.8I) . 

We will limit ourselves to considering the first and second moments of 
two quantities that are frequently encountered in applications and that have 
therefore been the subject of much earlier work, viz. the cell perimeter P 
and the cell area A. These are explicitly given by 

P = i? av (4A)- 1/2 A, A = Rl(AX)- 1 F 2 , (2.10) 

with the angular factors 

1 n 

F k = -J2p k m(^-? m + t&nf3 m+1 ), k — 1,2. (2.11) 

m=l 

Setting successively F = P, P 2 , A, A 2 we find that the corresponding inser- 
tions in the numerator of Eq. (12.61) are 

i P = [r^ + D/rH]^- 1 / 2 ^, 

I p2 = nW^F 2 , 
I A = \nW~ l F 2 , 

I A2 = \ n (n + 1)W~ 2 F 2 . (2.12) 
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The simulation steps for finding the numerator of Eq. (12.61) are the same as 
for p n except that (iv) and (v) are replaced with (iv') and (v') given below. 

(iv') We multiply the insertion Ip of the quantity F of interest by the 
weight exp(— V) and accumulate the value thus obtained. 

(V) In the end the total accumulated value is divided by the total number 
of attempts and by the estimate obtained for p n . This provides an estimate 
for (F) n . The numerical data shown will all be for A = 1. 

3 Results and discussion 

3.1 The distribution of V and the indicator 

Before discussing our results for the sidedness probability p n we briefly con- 
sider the quantities V and 6 that via (11.31) and (11.41) enter into its definition. 
Let -P(V) denote the probability distribution of V and <p n = (9) the proba- 
bility for an attempted cell generation to be successful. In terms of these we 
may rewrite (11.41) as 



which exhibits the important intermediate role of P(V) and n . 

In order to show what P(V) looks like, we have plotted its logarithm 
in Fig. [2] for n = 50, 100, 200, and 400. The curves clearly demonstrate 
that for n —>■ oo there is convergence to a limit. For V — > ±oo the limit 
distribution appears to decay exponentially, P(V) ~ exp(— k±|V|), but with 
very different decay constants: we obtain k + = 0.185 ± 0.05 from a fit in the 
range 3 < V < 30 followed by extrapolation ton = oo, and k_ = 2.47 ± 0.02 
from a fit in the range —3 < V < —1.5. This large V behavior has not yet 
been explained theoretically. 

The G function in (11.41) imposes constraint (12. 3p and is at the origin of 
the failed generation attempts. Whereas these do not contribute to p n in 
step (iv) of the algorithm above, Eq. (13.11) shows that via (6) = <p n they 
do enter into the determination of its normalization. In the last column of 
Table [1] we list the fractions (f) n of successful attempts as determined from 
the simulation. Although <ft n is equal only to 03 = 0.058 for n — 3, it turns 
out to rise rapidly with n, is already as high as pio = 0.8 for n = 10, and 
tends to unity for n — > oo. That is, attrition disappears in the large n limit. 

This brings out the two key steps that are responsible for the success of 
the present algorithm: (i) the limit distribution of P(¥) has become n inde- 
pendent since we extracted from the initial expression for p n the appropriate 
n dependent prefactor p^ given in (jl.2p ; and (ii) attrition disappears for 




(3.1) 



8 








Figure 2: Logarithm of the probability distribution P(V) of V [see Eq. (|3.ip ] for 
four different values of n, showing convergence to a limit distribution for n = oo. 

large n because of our choice of the angles (£, rf) as the variables of integra- 
tion. 

3.2 Sidedness probability p n 

In Table [1] we present the results for the sidedness probability p n for n in 
the range between n = 3 and n = 1600. They are based on a number N n of 
generation attempts given in the fourth column of that table. 

The second column of Table [1] shows the best results for p n found in 
the literature for each value of n. The p% value was obtained by numerical 
integration, the other p n are Monte Carlo results. For p 4 , . . . , p 7 the statistical 
error is in the last decimal; for higher n the standard deviations are indicated. 

Our own results for p n , given in the third column of Tabled! are accurate 
up to absolute errors of order less than 10~ 5 . Standard deviations were cal- 
culated by subdividing the data into twenty or more groups and considering 
the dispersion of the group averages. We will now discuss these results as a 
function of n. 
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Case n = 3. - The probability p^ for a cell to be three-sided is the only 
one that has been evaluated by numerical integration. This was done by 
Hayen and Quine [6], who reduced the original integral to a four-dimensional 
one. They present a 12-digit result of which the second column of Table [T] 
shows only the first seven significant decimals. For n = 3 we performed an 
especially long run with the purpose of testing our Monte Carlo method and 
checking the result of Ref. [6j. As shown in Tabled], our method reproduces 
six significant digits of Hayen and Quine's result and leaves their value within 
our error bars. 

For all n > 3 the literature results are based on Monte Carlo evaluation. 

Cases n = 4 through n = 7. The most accurate literature results in 
this intermediate regime are due to Calka [TJ, whose algorithm like ours fixed 
n in advance. Our results are fully compatible with those of Ref. [7J. 

Cases n = 8 through n = 15. Monte Carlo results obtained in the 
1980's by Brakke [8 J for 3 < n < 16 stayed for a long while the most accurate 
ones that were available. Our simulations confirm all of Brakke's results. 
Beyond n « 10 the accuracy of the Monte Carlo algorithm of Ref. [8 J rapidly 
goes down with increasing n, and for n — 16 its error bars are as large as the 
result itself. This effect is simply due to the low relative frequency of cells 
of so many sides, the number n not being fixed in this method. By contrast, 
the accuracy of our method, for a fixed amount of computer time invested 
per value of n, stays roughly constant. 

Case n > 16. Drouffe and Itzykson [9] developed a more powerful 
simulation method aimed at simulating cells of larger sidedness. In their 
method n is again fixed in advance. Their accuracy amounts to roughly a 
single significant digit in the regime 16 < n<25; forn>25 the error becomes 
again of the order of p n itself. This error increase is due to attrition, i.e., 
an increasing rejection rate of configurations that are generated but do not 
satisfy the required geometrical constraints. From our data it appears that 
for n ^ 25 the work of Ref. f9] misses the true values by an ever larger factor 
and that only their logarithmic order of magnitude is right |17j . Again, the 
method of the present work maintains an error in the fifth digit, i.e., a relative 
error not larger than 10~ 4 . 

Case of extremely large n. - The range of n from 50 to 1600 had so far 
been unexplored territory. In this very large n regime attrition is negligible 
and, for a constant calculational effort per value of n, the method keeps 
producing results with an error only in the fourth significant digit. The 
probabilities p n are extremely small. Numerically we could easily handle 
such small numbers by first factorizing out pn^ of which we computed and 
stored only the logarithm. As discussed in subsection 13.11 the remaining 
factor C n = (Oe _v ) has a finite distribution and hence causes no underflow 
problems. 

Generating the values of such 'unphysically' small probabilities is much 
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more than a mere technical achievement. First, it provides another check 
that the program works correctly; indeed we find that for n — > oo the ratio 
p n /Cpn = C n /C tends to unity as it should. Second, it gives access to the 
large n expansion of p n to be discussed in subsection 13.31 Thirdly and most 
importantly, values of n this large are required to see the separation of length 
scales that occurs in the many-sided cell; this is the subject of section HI 

Sum rules. - The p n should obey certain sum rules. Upon summing the 
p n of Table [1] and writing X n = Yln°=3 X n p n we find 

oo 

^Pn = 1.000 010(15), 

n=3 

n = 6.0001(1), 

r? = 37.7816(7), 

fi = ^-n 2 = 1.780 4, (3.2) 

with an error in fi at most equal to ±0.0015 but probably smaller due to par- 
tial cancellation of the errors in n 2 and n 2 . The first and second relations of 
(13.21) may be compared to the exactly known values 1 and 6, respectively. The 
second moment n 2 has an exact expression as a double integral [IB] , which 
when evaluated numerically gives n 2 = 37.780 811.... Hence [i = 1.780 811... 
numerically exactly. We therefore see that, when their error bars are taken 
into account, our Monte Carlo data are in excellent agreement with these 
sum rules. 

Conclusion. - The general conclusion of this subsection is that for low 
n (say n ^ 8) our method is probably as good as several of the existing ones. 
If we did slightly outperform them in that small n regime, that was only 
because of the length of our runs. However, for larger N (say n ^ 8) our 
method has a decisive advantage over the existing ones. 



3.3 Asymptotic behavior of p n 

On the basis of the numerical data we will now discuss the asymptotic be- 
havior of p n for n — > oo. Analytically it is known that p n = CnPn^ with p^ 
given by Eq. fll.2p and where the correction factor C n may be obtained by 
a series expansion that classifies contributions according to their power in 
n~ l l 2 . On that basis Ref. [3] fitted the limited p n data available at that time 
(essentially n ^ 30) by a correction term proportional to n^ 1 / 2 . It remained 
possible, however, that the coefficient of the n~ 1 / 2 term would cancel, and 
indeed Drouffe and Itzykson [H] had hypothesized earlier that the leading 
correction was of order ri~ l . 

The numerical data of this work now indicate unambiguously that the 



11 



series is actually one in powers of n 

C (8tt 2 ) ? 
Pn ~ 4vr 2 (2n)! 

where e%, €2, ■ ■ ■ , are numerical coefficients. 

The factor in square brackets in Eq. (I3.3P is equal to C n /C, which for 
n — > 00 is known to tend to unity. In Fig.[3j in order to find the corrections to 
the leading order term in (13. 3p . we have plotted n{l — C n /C) = ei—e 2 n~ 1 + . . . 
against n~ l . This figure shows that the intercept with the vertical axis is 
located at e\ = 14.00 ± 0.05. We may now proceed by subtracting this 
estimated value of e% from the curve of Fig. [3], multiply it again by n, and 
look for a new intercept with the vertical axis which, if it is well-defined, 
is equal to — e 2 . Upon iterating until the statistical errors obscure a well- 
defined intercept we obtained in this way estimates for the first few e^. The 
uncertainties increase with the index i. We found 

e x = 14.00 ±0.05, e 2 = 94 ± 4, e 3 = 375 ± 80, (3.4) 

in which the errors are correlated: the values deviate together upward or 
downward. The important conclusion is that p n has a series expansion in 
powers of n~ x . The cancellation of the half-integer powers in the expansion 
of Ref. pQ is no doubt due to a symmetry in the theory that still remains to 
be identified. 

In a final remark we wish to stress that finding this asymptotic expansion 
is different from finding a 'best fit', which we do not attempt here. The curve 
of Fig. [3] is close to the sum of a constant and an exponential in n' 1 , but we 
have no reason to believe that there exists a simple analytical expression that 
fits all data within their error bars. 

3.4 Perimeter and area 

We have simulated the two cell properties that have received the greatest 
attention in the literature, viz. the cell perimeter P and the cell area A. 
We determined the average, the second moment and the variance of both of 
these quantities as a function of n. The perimeter results are summarized in 
Table [2] and the area results in Table [3j 

Similar tables extracted from the literature were compiled by Okabe et al. 
[2J. However, by far the most accurate ones appear in unpublished work by 
Brakke [S] and concern the regime 3 < n < 16. All our area and perimeter 
data are compatible with those of Ref. [8], but our error bars are strongly 
reduced. A further check on the numerical data is provided by two more 
sum rules, 

7^ = 4.000 035(65), ~A~ n = 1.000 02(2), (3.5) 



1 _ £1 £2 _ £3 

00 



n 



n- 



(3.3) 
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n 


Literature 

Pn 


This work 

Pn 


N n 




3 


1.124001... x 10 


-2 


(1.124000 ±0.000021) x 10" 2 


1.2 x 10 10 


0.0580 


4 


1.06838 x 10 


-1 


(1.068454 ± 0.000025) x 10" 1 


2 x 10 9 


0.1730 


5 


2.5946 x 10 


-1 


(2.59444 ±0.00007) x 10" 1 


1.6 x 10 9 


0.3077 


6 


2.9473 x 10 


-1 


(2.94723 ±0.00009) x 10" 1 


2 x 10 9 


0.4391 


7 


1.9877 x 10 


-1 


(1.98768 ±0.00007) x 10" 1 


4 x 10 8 


0.5564 


8 


(9.0116 ± 0.0020) x 10 


-2 


(9.0131 ±0.0006) x 10" 2 


10 8 


0.6554 


9 


(2.9644 ±0.0012) x 10 


-2 


(2.9652 ±0.0002) x 10" 2 


8 x 10 7 


0.7361 


10 


(7.4471 ± 0.0059) x 10 


-3 


(7.4487 ±0.0006) x 10" 3 


8 x 10 7 


0.8002 


11 


(1.4796 ± 0.0026) x 10 


-3 


(1.4818 ±0.0002) x 10~ 3 


6 x 10 7 


0.8501 


12 


(2.409 ±0.011) x 10 


-4 


(2.4000 ±0.0002) x 10" 4 


6 x 10 7 


0.8884 


13 


(3.18 ±0.04) x 10 


-5 


(3.2324 ±0.0003) x 10" 5 


6 x 10 7 


0.9175 


14 


(3.60 ±0.13) x 10 


-6 


(3.6835 ± 0.0004) x 10" 6 


4 x 10 7 


0.9393 


15 


(3.7 ±0.4) x 10 


-7 


(3.6017 ±0.0004) x 10" 7 


4 x 10 7 


0.9556 


16 


(2.3 ±0.3) x 10 


-8 


(3.0574 ±0.0004) x 10" 8 


4 x 10 7 


0.9677 


17 






(2.2762 ±0.0002) x 10" 9 


4 x 10 7 


0.9765 


18 


(1.3 ±0.5) x 10" 


10 


(1.4989 ± 0.0002) x 10~ 10 


4 x 10 7 


0.9830 


19 






(8.7983 ±0.0013) x 10~ 12 


4 x 10 7 


0.9878 


20 


(1.5 ±0.8) x 10- 


13 


(4.6314 ±0.0004) x 10~ 13 


8 x 10 7 


0.9912 


21 






(2.1994 ±0.0004) x 10" 14 


2 x 10 7 


0.9937 


22 






(9.4835 ±0.0017) x 10" 16 


2 x 10 7 


0.9955 


23 






(3.7227 ± 0.0005) x 10" 17 


2 x 10 7 


0.9968 


24 






(1.3379 ±0.0003) x 10" 18 


2 x 10 7 


0.9977 


25 


(9.6 ±5.9) x 10- 


21 


(4.4184 ±0.0004) x 10" 20 


4 x 10 7 


0.9984 


30 


(1.3 ± 1.1) x 10" 


29 


(5.4595 ± 0.0005) x 10" 28 


4 x 10 7 


0.9997 


40 


2.4 x 10- 


50 


(6.7349 ±0.0006) x 10" 46 


8 x 10 7 


1.0000 


50 


1.5 x 10" 


75 


(5.223 ±0.001) x 10" 66 


1.6 x 10 7 


1.0000 


60 






(7.192 ±0.002) x 10" 88 


1.2 x 10 7 


1.0000 


70 






(3.4805 ± 0.0004) x 10" m 


3 x 10 7 


1.0000 


80 






(9.598 ± 0.002) x 10" 136 


10 7 


1.0000 


90 






(2.1616 ±0.0005) x 10" 161 


0.8 x 10 7 


1.0000 


100 






(5.2691 ±0.0006) x 10" 188 


1.6 x 10 7 


1.0000 


150 






(1.0535 ±0.0002) x 10" 332 


4 x 10 6 


1.0000 


200 






(3.818 ±0.001) x 10" 492 


4 x 10 6 


1.0000 


300 






(1.084 ±0.001) x 10" 841 


2 x 10 6 


1.0000 


400 






(9.863 ±0.003) x 10" 1221 


4 x 10 6 


1.0000 


600 






(3.645 ±0.002) x 10" 2040 


10 6 


1.0000 


800 






(1.326 ±0.001) x 10~ 2918 


2 x 10 6 


1.0000 


1000 






(6.365 ±0.003) x 10~ 3841 


1.6 x 10 6 


1.0000 


1200 






(3.262 ±0.002) x 10" 4798 


1.2 x 10 6 


1.0000 


1400 






(1.385 ±0.001) x 10" 5784 


0.8 x 10 6 


1.0000 


1600 






(7.4306 ± 0.0020) x 10" 6796 


4 x 10 6 


1.0000 



Table 1 : The sidedness probability p n . Second column: literature data taken from Hayen 
and Quine [6] for p^; from Calka [7] for p^, . . . ,pj; from Brakke [8] for p$, . . . ,pis; and 
from Drouffe and Itzykson [9] for p n with n > 16. Third column: p n and its standard 
deviation calculated by the Monte Carlo method of this work. Fourth column: number 
N n of cell generation attempts. Fifth column: fraction 0„ of successful attempts. 
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n 


(P)n 




(P 2 )n 




(SP 2 ) 1 / 2 


3 


2.740296 


(2) 


8.17130 


(2) 


0.81368 


4 


3.219524 


(3) 


11.04819 


(2) 


0.82634 


5 


3.642658 


(3) 


13.96626 


(3) 


0.83504 


6 


4.026307 


(4) 


16.91958 


(4) 


0.84169 


7 


4.380000 


(6) 


19.90272 


(6) 


0.84702 


8 


4.710196 


(8) 


22.91084 


(8) 


0.85140 


9 


5.020869 


(12) 


25.94026 


(12) 


0.85506 


10 


5.315211 


(11) 


28.98790 


(12) 


0.85816 


11 


5.595488 


(10) 


32.05043 


(12) 


0.86078 


12 


5.863536 


(11) 


35.12588 


(13) 


0.86304 


13 


6.12072 


(2) 


38.2114 


(2) 


0.86497 


14 


6.36824 


(1) 


41.3055 


(2) 


0.86664 


15 


6.60705 


(2) 


44.4066 


(2) 


0.86809 


16 


6.83797 


(2) 


47.5136 


(3) 


0.86936 


17 


7.06173 


(2) 


50.6258 


(2) 


0.87047 


18 


7.27884 


(2) 


53.7410 


(3) 


0.87145 


19 


7.48992 


(2) 


56.8598 


(3) 


0.87233 


20 


7.69544 


(1) 


59.9820 


(3) 


0.87310 


21 


7.89576 


(3) 


63.1066 


(4) 


0.87380 


22 


8.09118 


(3) 


66.2318 


(5) 


0.87442 


23 


8.28215 


(2) 


69.3596 


(3) 


0.87499 


24 


8.46892 


(3) 


72.4890 


(4) 


0.87551 


25 


8.65171 


(2) 


75.6198 


(4) 


0.87598 


30 


9.51379 


(2) 


91.2825 


(5) 


0.87783 


40 


11.03971 


(1) 


122.6501 


(4) 


0.88005 


50 


12.37983 


(2) 


154.0370 


(5) 


0.88135 


60 


13.58887 


(3) 


185.4355 


(7) 


0.88219 


70 


14.69896 


(2) 


216.8384 


(6) 


0.88279 


80 


15.73105 


(2) 


248.2460 


(7) 


0.88323 


90 


16.69951 


(3) 


279.6545 


(8) 


0.88357 


100 


17.61487 


(2) 


311.0645 


(8) 


0.88384 


150 


21.61817 


(3) 


468.1280 


(12) 


0.88465 


200 


24.98833 


(3) 


625.1998 


(11) 


0.88505 


300 


30.63607 


(2) 


939.3527 


(12) 


0.88544 


400 


35.39384 


(2) 


1253.5092 


(13) 


0.88564 


600 


43.37090 


(4) 


1881.820 


(3) 


0.88584 


800 


50.09346 


(2) 


2510.143 


(2) 


0.88593 


1000 


56.01490 


(1) 


3138.456 


(2) 


0.88599 


1200 


61.36764 


(1) 


3766.773 


(2) 


0.88603 


1400 


66.28953 


(2) 


4395.087 


(4) 


0.88606 


1600 


70.87047 


(2) 


5023.408 


(2) 


0.88608 


oo 










0.886226... 



Table 2: Estimates of the average (P) n , the second moment {P 2 ) n , and the root-mean- 
square fluctuation (5P 2 )l/ 2 of the cell perimeter P. The numbers in parentheses represent 
the standard deviation in the last digit. The entries of the third column have an error of 
at most one unit in their last digit. The limit value ^tt 1 ^ 2 — 0.886226... for n = oo has 
the status of a conjecture. 
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n 


(A)n 




(A 2 )n 






3 


0.343087 


(3) 


0.161573 


(3) 


0.12092 


4 


0.558052 


(4) 


0.401285 


(5) 


0.14989 


5 


0.774080 


(4) 


0.73675 


(!) 


0.16586 


6 


0.995789 


(5) 


1.17953 


(2) 


0.17698 


7 


1.22251 


(!) 


1.73516 


(3) 


0.18541 


8 


1.45328 


(!) 


2.40724 


(3) 


0.19200 


9 


1.68736 


(!) 


3.19847 


(4) 


0.19756 


10 


1.92408 


(2) 


4.11064 


(7) 


0.20213 


11 


2.16295 


(2) 


5.1451 


(!) 


0.20599 


12 


2.40366 


(2) 


6.3033 


(!) 


0.20929 


13 


2.64578 


(2) 


7.5854 


(!) 


0.21217 


14 


2.88906 


(3) 


8.9920 


(!) 


0.21469 


15 


3.13331 


(3) 


10.5234 


(!) 


0.21689 


16 


3.37835 


(3) 


12.1797 


(2) 


0.21885 


17 


3.62416 


(2) 


13.9619 


(2) 


0.22059 


18 


3.87034 


(3) 


15.8680 


(2) 


0.22215 


19 


4.11703 


(3) 


17.8996 


(3) 


0.22357 


20 


4.36415 


(3) 


20.0570 


(3) 


0.22484 


21 


4.61158 


(4) 


22.3394 


(3) 


0.22601 


22 


4.85923 


(4) 


24.7464 


(4) 


0.22706 


23 


5.10715 


(4) 


27.2790 


(4) 


0.22803 


24 


5.35531 


(5) 


29.9371 


(5) 


0.22891 


25 


5.60358 


(6) 


32.7198 


(6) 


0.22974 


30 


6.84686 


(6) 


48.5090 


(6) 


0.23306 


40 


9.33913 


(4) 


89.470 


(!) 


0.23723 


50 


11.83458 


(7) 


142.932 


(2) 


0.23977 


60 


14.33183 


(6) 


208.900 


(2) 


0.24146 


70 


16.82979 


(6) 


287.363 


(3) 


0.24267 


80 


19.3283 


(!) 


378.331 


(4) 


0.24358 


90 


21.82726 


(8) 


481.800 


(4) 


0.24429 


100 


24.32627 


(6) 


597.763 


(3) 


0.24485 


150 


36.8236 


(2) 


1365.10 


(!) 


0.24658 


200 


49.3224 


(!) 


2444.94 


(!) 


0.24742 


300 


74.3214 


(2) 


5542.16 


(3) 


0.24826 


400 


99.3206 


(!) 


9889.33 


(3) 


0.24871 


600 


149.3196 


(3) 


22333.58 


(!) 


0.24914 


800 


199.3201 


2) 


39778.23 


(6) 


0.24935 


1000 


249.3192 


(2) 


62222.3 


(!) 


0.24948 


1200 


299.3193 


(2) 


89666.7 


(!) 


0.24957 


1400 


349.3187 


(3) 


122110.8 


(2) 


0.24963 


1600 


399.3190 


(2) 


159555.4 


(2) 


0.24968 


00 










0.250000... 



Table 3: Estimates of the average (A) n , the second moment {A 2 ) n , and the normalized 
root-mcan-square fluctuation n~ 1//2 (SA 2 }1/ 2 of the cell area A. The numbers in parentheses 
represent the standard deviation in the last digit. The entries of the third column have an 
error of at most one unit in their last digit. The limit value \ for n = oo has the status of 
a conjecture. 
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0.1 . . 0.2 0.3 



Figure 3: To study the asymptotic large n behavior of the sidedness probability 
Pn = C n pn ^ we plot the quantity n(\ — C n /C) where C = lirrin^oo C n [see Eqs. (01) 
and p.3p ]. The solid line connects the data points. Data shown are in the range 
3 < n < 300. The largest error bars occur for small 1/n and are of the order of the 
data symbols. The intercept of the curve with the vertical axis is the coefficient 
ei of the leading correction term in the expansion of Eq. f)3.3() . 
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0.01 0.02 

l/n 



Figure 4: The asymptotic large-n behavior of the root mean square fluctuations 
(SP 2 )l/ 2 and (SA 2 )l/ 2 of the perimeter and the area, respectively, of the n sided 
cell. Data shown are in the range 50 < n < 1600. The error bars are smaller than 
the data points. The two curves that connect the data points are asymptotically 
straight lines which for n — > oo both appear to converge ^. 
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for which the exact values are 4 and 1, respectively. 

We now turn to the large n behavior. Our data indicate the expansions 

/ \i -i 

(imjz+ain 2 + ci3 n 2 + ..., 

2 2 

nn + b + bin^ 1 + . . . , 
\n + c + ci^" 1 + . . . , 

(±n) 2 + cL 1 n + d + ... , (3.6) 

which again go down by integer powers of n. They imply that the variances 
have the series 

(P 2 )n-(P)l = b -2^a h + ... 

(A 2 ) n -(A) 2 n = (d_ 1 -I Co ) n + ... . (3.7) 

The leading terms in each of the four series of Eq. (13. 6 h are known from 
theoretical analysis [3j [1] . Heuristically they follow from the sole observation 
that in the large n limit the Voronoi cell becomes a circle of radius R c = 
in/An) 1 ! 2 . Theoretical analysis can in principle also produce the higher order 
terms in (13. 6p . but this has not been attempted yet. Consequently, the 
leading coefficients of the series in (13. 7p are not known analytically. Here our 
simulations results provide answers. 

In Fig.Hwe have plotted [(6P 2 ) n /7i] 1/2 and 2[(5A 2 ) n /n] 1/2 . The numerical 
data strongly point to limit values equal to | for both quantities when n — > 
oo. Conjecturing that these limits are exact we then conclude that 

b - 2n^ai = |tt, d_i - |c = ^. (3.8) 

Analysis of the (P 2 ) n data to next order suggests that (P 2 ) n / (nn) — 1 tends 
to — 1 as n —>■ oo. Conjecturing that this, too, is exact and combining it with 
the first one of Eqs. (13. 8 ft we arrive at 

ai = -§7r*, b = -7r, (3.9) 

We do not attempt, however, a similar analytical conjecture for the second 
pair of coefficients, Co and nor will we pursue estimates for the higher 
order coefficients in the series (13.61) and (13.71) . except below in connection 
with Lewis' law. 

Lewis' celebrated law [19j is an empirical statement about one of the 
cell's most conspicuous properties, viz. the relation between its area and its 
number of sides. The law states that the average area (A) n of an n-sided cell 
increases with n as 

(A) n =^(n-n ), (3.10) 



(P)n = 

(P\ = 

(A)n = 

(A\ = 
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where ao and uq are constants and where we have displayed again the de- 
pendence on the areal particle density A. Sometimes (see the discussion in 
Ref. pj) this law is written in the more restricted one-parameter form 

(A) n = - A 6) - 1 (3.11) 

It is found, however, that (A) n deviates from linearity with n in simulations 
of Poisson-Voronoi diagrams as well as in the data from most experimental 
systems. We now look at what the asymptotic analysis has to say 
In Refs. [31 H] we proved that asymptotically 

(A) n ~7rR 2 c = —, n^oo. (3.12) 

and this result has been incorporated in the series for (A) n in (13.61) . A 
coefficient ao ~ \ had since long been suspected by various authors [HI [201 12] • 
Going now beyond (I3.12p and proceeding in the same way as for p n , we can 
determine the coefficients of the series of (13.61) for (A) n on the basis of our 
simulation results of Table [2j This yields 

Co = -0.6815(5), d = 0.750(5), c 2 = 3.15(10). (3.13) 

We now consider the laws (I3.10p and (13. lip . The fact that we found ci, C2, . . . 
to be nonvanishing confirms once more that (A) n is not strictly linear in n. 
From the above it follows that in (I3.10p one should choose 

n = _4c = 2.7260(4) (3.14) 

if one wishes it to correctly represent the asymptotic behavior of {A) n for 
Poisson-Voronoi cells. This is of course different from finding a best fit to a 
limited set of (A) n data in a restricted n interval. If that is the purpose, the 
values of ao and no will depend on the available data and on the way the fit 
is carried out. The one-parameter law (13.111) postulates a relation between 
ao and no that is violated in the asymptotic expansion. Hence (13.111) cannot 
be used to describe the large n behavior of (A) n and merely has the status 
of an empirical fit to the data, the value of b again depending on the data 
set and on how the fit is done. 

4 Characteristic cell shapes 

It has been established [H [1] that in the large n limit the Voronoi cell tends 
to a disk of radius R c = (n/A-n) 1 / 2 [22]. In Ref. pQ it was furthermore shown 
that the cell perimeter undergoes 'elastic' deformations from circularity, the 
elasticity being, of course, of entropic origin. The probability law of these 
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deformations was given in the large n limit. In this section we show how 
our Monte Carlo method allows us what was hitherto impossible, namely to 
simulate for any finite n the detailed statistics of the cell shape. 

We Monte Carlo generated cells of prescribed sidedness n in a 'natural', 
that is, an unbiased, environment. This was done as follows. For a given n 
the cell angles (£, rf) were drawn randomly and /3* was found according to 
the rules of section 12 .31 The cell radius was taken equal to its most probable 
value R c = (n/47r) 1 / 2 and the cell boundary was constructed. This boundary, 
together with the position of the central particle, fixes the positions of the n 
first neighbor particles. We then determined the cell's fundamental domain 
J 7 , that is, the union of the n disks of radius S m centered at the vertices S m . 
The complement of T in a large rectangle of suitable size was subsequently 
filled randomly with particles of a uniform density equal to one. The particles 
added by this procedure are all necessarily second or higher order neighbors 
of the central one. The Voronoi construction was finally applied to the full 
collection of particles to complete the Voronoi cell diagram. 

4.1 Cells of n = 24, 48, and 96 sides 

We have generated typical cell shapes for a sequence of values of n, starting 
with n = 3 and doubling n each time. Figs. El EJ and [7J in which the dots 
represent the particles, show the results for cells of n — 24, 48, and 96 
neighbors. The three pictures are at different scales but all have unit particle 
density. This picture sequence illustrates the tendencies that characterize 
many-sided cells. One tendency is for the first neighbor cells to be elongated. 
This feature is apparent already for n = 24 and becomes very pronounced for 
n = 48, whereas the n = 96 cell has only very elongated neighbors. The same 
phenomenon was observed by Lauritsen et al. (21] , but in a different system. 
These authors consider Poisson- Voronoi diagrams to which they assign an 
'energy' that favors many-sided cells. Snapshots of their configurations show 
a dense structure of many-sided cells (of sidednesses higher than n = 60) 
separated by mostly four-sided elongated cells. Their procedure does not, 
however, provide estimates for p n in an unbiased Poisson- Voronoi diagram. 

Another tendency, appearing similarly in Ref . [21] , is for the first-neighbor 
particles to align on what tends towards a continuous curve. Whereas for 
n = 24 some imagination is still needed to see this curve, it becomes clearly 
distinguishable for n = 48 and is immediately obvious to the eye for n = 
96. The typical distance between nearest neighbor particles along this curve 
decreases as 2tt(2R c ) /n ~ n^ 1 / 2 . We note that whereas Voronoi cells are 
always convex, the 'curve' connecting the first neighbors need not enclose 
a convex area; in each of the Figs. El El and [7] there are small but clearly 
distinguishable deviations from convexity. 

Fig. [8] enlarges a detail of Fig.[7] and shows a collection of first-neighbor 
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cells. All first neighbors fully visible in the figure are four-sided except those 
marked A, B, D, E, which are five-sided, and the cell C, which is either five- 
or six-sided (this depends on how the two almost coinciding three-vertices 
are arranged at the point marked ! 2V; a higher resolution is needed to decide 
this question). Fig. [H] illustrates that in the limit of large n four-sided first 
neighbors become dominant. In Ref. [I] it was argued that five-sided cells 
constitute a fraction only of order n^ 1 ^ 2 of all first-neighbor cells, and that 
the probability of six- and higher-sided first-neighbors is of still higher order 
in n _1//2 . In Fig.[8]the cell marked P is a second neighbor to the central cell. 
The boundary separating it from the first neighbors has been shown as a 
heavy solid line on which we will further comment shortly. 

4.2 Very large cells 

Focusing now on very large n we show in Fig.[9]a central particle located in the 
origin and having 1536 neighbors. As before, the dots represent the particles. 
The inner contour, which is nearly indistinguishable from a circle of radius 
R c , represents the boundary of the Voronoi cell of the central particle. The 
outer 'curve', which is also very close to circular but of radius 2R C , represents 
the alignment of the 1536 first-neigbor particles. Their high line density gives 
the impression of a continuous curve. Cell boundaries other than those of 
the central cell have not been drawn; they would totally blacken the empty 
annular region between the boundary of the central cell and its first-neighbor 
particles. 

The boxed region in Figgis shown enlarged in Fig.[TUJ where we did draw 
all Voronoi cell boundaries. The extreme elongation of the first-neighbor 
cells is what first strikes the eye. The discrete structure of the 'curve' of 
first neighbors is also apparent now. The distances t m between successive 
first-neighbor particles along this curve are of order n _1//2 . More precisely, if 
we set i m = A m (47r/n) 1//2 , then the theory pQ implies that for n — > oo the X m 
are independent identically distributed random variables of probability law 
A m exp(— A m ). Random deviations from a local straight line are too small 
to be discernible to the eye; they may be argued [26] to decrease as n~ 3 / 2 , 
which is also the order of magnitude of the systematic deviations due to the 
radius of curvature 2R C . The large cell and its environment are characterized, 
therefore, by four different length scales, each varying with its own power of 
n. They have been summarized in Table H] below. One has to go to n values 
as high as we did in order for the separation of scales to become clearly 
visible. 

Very large n is required also for still another feature to become apparent. 
In Fig.[Tn]the boundary between a second-neighbor cell and its adjacent first- 
neighbor cells is, by construction, composed of points equidistant to the 
second-neighbor particle and the almost continuous straight line of first— 
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Scale Length 

n 1 / 2 Cell radius 

n° Typical interparticle distance outside the first neighbor circle 

?i -1 / 2 Typical distance between successive first neighbor particles 

n -3/2 Random deviations of first neighbors from full alignment 

Table 4: Four length scales characterizing the n-sided Voronoi cell in the large n limit; 
see Figs.HandrjDl 



neighbor particles. But such a boundary is a parabola. Hence, in the limit 
n — > oo the boundary separating the set of first from the set of second-neigbor 
cells is piecewise parabolic. Indeed, with this observation in mind one now 
recognizes the boundary segment of cell P in Fig.[H] (heavy solid line), and 
others in that same figure, as 'incipient parabolic'. Such knowledge was at 
the basis of the theory of two-cell correlations exposed in Ref. @]. There, 
laws discovered in the n —>■ oo limit were extrapolated backward and shown 
to be relevant for finite n. It was shown, in particular, that Aboav's linear 
relationship [121 [2] between n and the total average sidedness nm n of the 
neighbors of an n-sided cell cannot hold in Poisson- Voronoi diagrams. We 
expect that in the future the study of large cells will shed further light also 
on various issues relevant for the finite n behavior. 

5 Summary and conclusion 

In this paper we have developed a new Monte Carlo method for evaluating the 
sidedness probability p n for arbitrary n. The method, which is constructed 
on the basis of an extension of earlier theory [3J, [1] , is not difficult to imple- 
ment once one has available the rather complicated analytic expressions that 
intervene. 

We have determined p n as well as the first and second moments of the 
n dependent cell perimeter and cell area. Full agreement is obtained with 
earlier results for p n due to Hayen and Quine |6J, Calka [7J, and Brakke [8], 
whose data extend up to n = 16. For n^lO we have reduced the error bars on 
p n very considerably. In the range up to n = 50 we improved and corrected 
the p n data due to Drouffe and Itzykson [9]. For 50 < n < 1600 we obtained 
data in a range that that had so far remained inaccessible. This enabled us 
to investigated the asymptotic large n behavior of p n and of the perimeter 
and area moments. On the basis of our numerical results we conclude that 
they all have asymptotic series in entire powers of n -1 , possibly up to an 
overall pref actor na. 

Exploiting our full control of the cell statistics we have exhibited occur- 
rences of extremely rare many-sided cells in a typical environment of ordinary 
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-2.5 2.5 

Figure 5: A typical Voronoi cell with n = 24 neighbors. The dots represent the 
particles. 
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24 




25 



Figure 8: Enlargement of the box in Fig. [7] showing some of the strongly elongated 
first neighbors of the central cell. Among the first neigbors fully visible, cells 
A, B,C, D, and E have more than four sides. The arrow marked '2V points 
to two three-vertices that coincide within the resolution of the figure. Cell P is 
a second neighbor whose boundary with the first neighbors (heavy lines) is an 
example of an 'incipient parabola segment' (see text). 
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Figure 9: Approach to the infinite n limit. The origin is occupied by a particle 
whose Voronoi cell has n = 1536 sides. The almost circular inner curve is the cell 
boundary of the central cell. The other cells boundaries have not been drawn. The 
almost circular outer curve is made up of the 1536 first-neighbor particles. The 
region inside the box is shown enlarged in Fig.llOl 
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Figure 10: Enlargement of the box in Fig.0 where now all cell boundaries have 
been drawn. The discrete structure of the outer 'curve' of Fig.[9]has become visible 
here. 
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cells. Their embedding involves four distinct length scales, varying with four 
different powers of n. This has confirmed, among several other things, the 
very elongated shape of the first-neighbor cells. 

No particular effort was made to optimize the code. Our total investment 
of computer time on a recent model PC was limited to a few hundred hours 
and allowed us to obtain p n to at least four or five significant decimals for the 
set of n values listed in the tables. We have also not attempted to provide 
any 'best fits' to the numerical curves, as we have no reason to believe that 
there exist simple analytic expressions that fit the data within our error bars 
over their full range. 

The Monte Carlo work of this paper became possible only after initial 
analytic progress [31 [1] . We believe that it will in return spur further analytic 
investigation. One branch of study may concern the nature of the asymptotic 
expansions uncovered here. Another one may deal with correlations between 
a cell and its second, third, and higher topological neighbors, which are a 
recurrent theme in the theory and applications of Voronoi cells. 
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A Theory 

We present here the extension of earlier work that opens the way to the 
numerical simulations of this work. We consider uniformly and independently 
distributed particles in the plane. Let a particle be placed in the origin and 
let n other particles occupy the positions 2R 1; 2R 2 , . . . , 2R n . The sidedness 
probability p n of the cell containing the origin may then be written as a 
2n-fold integral on the midpoint coordinates [9j [TOj [131 [7], 

p n = 1 f dRi . . . dR n . . . , Rn) e -^ Rl -- R "). (A.l) 

n\ J 

Here the indicator function x is equal to unity (or to zero) on the domain 
of phase space where the perpendicular bisectors of the vectors 2R m , for 
m = 1,2, ... ,n, define an n-sided (or a fewer-sided) cell around the origin; 
and A is the two-dimensional volume of the area that should be void of 
particles if this cell is not to be intersected by any of the bisectors of the 
position vectors of the remaining particles. Explicit expressions for A and x 
are given in Refs. [3 [TJ. 
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A.l Starting point 

After one integrates over a common radial scale, expression (lA.lj) takes the 
form [Tj of an integral on the angle fli and on the sets of angles £ = {£ m } 
and 7] = {r] m }, 

Pn = / dft /d£dr/ 8( V £ TO - 2tt) <f( V Vm ~ 2tt) 

271 J ~*/ 2 J m=l m=l 

[jJ^ Tm ] [7r(1 + n -V )r , (A.2) 



m=l 



where G and /3* are as defined in the main text [Eqs. (12.21) and (12.11) ] and the 
derivative G' = dG/dfti is given explicitly by (IB. 81) and (IB. 41) of Appendix 
IB1 The definitions of the new symbols occurring in (1A.2|) follow below. 



First of all, the symbol J d£dr/ in (IA.2I) is shorthand for the nested inte- 
grations 



d£d V = d6 / d Vl / d£ 2 ... 

Jo Jo Jo 

... / dry n _i / d£ n / drj n . (A. 3) 

Jo Jo Jo 



The notation is hybrid; the variables 71, /?2, 73, . . . , /3 n occurring here should 
be viewed as functions of £, r/, and the 'angle of rotation' f3\. They are given 
by 



m— 1 



£=1 
m— 1 

7m = -P*(€,V) + -ffe) + 771=1,..., 71, (A.4) 



where Y^e=i denotes an empty sum. Next, the T m and p m are functions of 
the 7 m and (3 m given by 

r m = sm ^ m \ + n lm \ m = l,...,n, (A.5) 
cos 2 (3 m 

cos7 m cos7 m _i . . .COS71 

Pm = n n Pn , 771 = 1, . . . , 71 - 1, (A.6) 

cos p m cos p m _i . . . cos pi 
and the condition 

n 

- iy £pm=±- (A.7) 



71 

m=l 
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Finally V is given by 

n 

V = — 2j P 2 m [ tan 7m - 7m + tan /3 m+ i - f3 m+ i 

m=l 

+ 7 m tan 2 7 m + /3 m+1 tan 2 /3 m+ J 

n 

+^I>m- 1 )(7m+/w)- ( A - 8 ) 



m=l 



The factor n included in its definition makes that, typically, V is of order n° 
as n —>■ 00. This completes the definition of the multiple integral flA.2p for 

Pn- 



A. 2 Transformations 



We now depart from the development of Ref. [I] and transform expression 
flA.2j) as follows. We integrate over f3\ and henceforth when writing (3\ it will 
be understood that it takes the value (3\ = /3* The integration requires 

that Eq. (12. ip have a solution. In Ref. p] a unique solution was shown to exist 
perturbatively for large n; in Appendix B of the present work we provide the 
demonstration for general n. We furthermore replace the upper integration 
limits of the integrals over the £ m and r\ m by 00 at the expense of introducing 
Heaviside theta functions. Using that £ m — f3 m = 7 m and r\ m — 7 m = p m +i 
and introducing a factor #(| — Pi), which may be done for free we find that 
Eq. (IA.20 may be converted into 



(n-1)! 
2mr n 



df 1 f 1 • • • df n f n / dr}i...dr], 



x 



5(J2 ^ ~ 2tt) 5(J2 Vm ~ 2n) Qe 



m=l 



m=l 



(A.9) 



in which 



and 



n 

G'{^ m P*)- 1 [ II Pm^mC 1 ] (1 + n~ l V)- n (A.10) 



m=l 



m=l m=l 

Expression ( 1A.9j) is more symmetric than ( 1A.2l) -( TA~3l) . Its integrand is a 
function exclusively of the £ m and ?] m . We have purposefully included extra 
weights £ m in the integrations in (1A.9I) and compensated for these by factors 
in the product on m in (1A.10I) . In this way we obtain property that T m £ m 
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remains finite when £ m — > 0, which was desirable analytically [3l [T] and is 
also necessary numerically. 

The same quantity V as defined in (1A.10I) was studied analytically in 
Refs. [21 [I] , where it was shown that for n — > oo it remains, typically, of order 



n°. 



(A.12) 



One further rewriting is useful. We set 

£m = «2m-l + a 2m, m = l,2, . . . , 71, 

and use that 

/•oo /*oo 

/ da2m-lda 2 m /(«2m-l + «2m) = / d£ m £mf(£m) (A.13) 

for any function /(£ m ). This converts (1A.9I) into the final result 

Pn = ^ 0) (Oe- v ), (A.14) 
where for any function X of the angular variables the average (X) is defined 

by 



(X) 



(°) , 

Pn -JO 



da i ... da 



2/i 



drjx... drj n 



2/i 



X 



5(J2 a m-27i)5(J2vm-27c)X. (A.15) 



m=l 



m=l 



,(0) 

!// - 1)! 



The normalization factor p n that appears here is easily calculated as 

2n 



2mr n 



'-Jo 



dai . . . da 2n S( « m - 27r) 



m=l 



d?7i . . . dr/ n S ( r/ m - 2tt) 



m=l 



(n-l)l (2tt)^- 1 (2tt 

X — — X 



in— 1 



2nvr n " (2ra-l)! (n-1)! 
87r 2 ) n 



4vr 2 (2n)! 



(A.16) 



which is (II .2p . This way of arriving at is slightly simpler than the original 
calculation of Ref. pQ. Expressions ( lA.14j) -( TA.i5l) are new and are at the basis 
of the Monte Carlo simulation of this work. The integrals in ( 1A.15I) directly 
suggest step (i) of the algorithm of subsection (12.31) . 
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B The equation G = 

We discuss here the function G defined by 

e 27rG _ COS71 cos 72 . . . COS 7 n 
COS Pi COS Pi . . . COS P n 



(B.l) 



The transformation to angular variables in Appendix El led to Eq. ( 1A.11I) . 
i.e. to the upper limits of integration P m -,lm < |. Since P m + 1m = £m and 
since £ m > 0, we have in fact that in the integral for p n the angles p m and 
7 m are restricted by 

- f < Pm, Ira < f ■ (B.2) 

Hence G — ► — oo whenever any of the 7 m tends to ±|, and G —>■ oo whenever 
any of the /3 m tends to ±|. We set now 

7m = 7m- A, m = l,...,n, (B.3) 

where the (3 m and 7 m are functions of the £ m and 7? m that may be read off 
by a comparison of Eqs. (IB. 31) and (IA.4[) . 



m—1 



1=1 

m—1 

7m = ^2(&-Vi)+£m, m = l,...,n. (B.4) 

^=i 

where again Yle=i denotes the empty sum. Making all Pi dependence explicit 
we get 

e 2nG = CQS(71 - Pi) CQS(7 2 - ft)... COs(% - fa) ^ 

cos(/3i + Pi) cos(/3 2 + . . . cos(/3 n + Pi) 

which we wish to study as a function of the single variable Pi, at fixed (£, 77). 
Expression (IB. 51) shows that exp(27rG) is positive on the interval 

- I + max 7 m < Pi < I — max /3 m , (B.6) 

l<m<n l<m<n 

provided this interval is not empty, that is, provided 

max 7 m + max /3 m < ir. (B.7) 

l<m<n l<m<n 

Because of the preceding discussion, G approaches —00 and 00 as Pi ap- 
proaches the left and right hand end points of this interval, respectively. To 
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show that G is actually monotonous in j3\ on the interval (1B.6j) . it suffices to 
analyze the derivative 



dG 



1 " 

- [tan(7 m - fa) + tan(/L + Pi)] ■ (B.8) 



m=l 

Since (3 m + 7 m = £ m > 0, it follows that there are three cases, namely (i) 
An, 7m > 0, (ii) p m > 0, 7 m < 0, and (iii) (3 m < 0, 7 m > 0. By considering 
each of them separately one deduces that the summand in Eq. (1B.8|) is always 
positive. It follows that dG/d/3i > and hence that G = has a unique 
solution Pi = /?*(£, ?y) in the interval ( IB. 61) . 

Hence we have shown that the conditions £ m ,r) m > and P m ,j m < § 
suffice for the equation G = to have a unique solution /?* in the physical 
interval (1B.6I) . This condition involves, however, the angles P m and 7 m which 
are determined by the solution We would like to have a criterion for the 
existence of a solution in terms of the sole sets (£, rj) given at the outset. 
By retracking the solution method, we see that it is valid for all (£, rj) as 
long as the 'physical' interval flB.6|) is not empty, that is, as long as Eq. (IB. 71) 



is satisfied. When made explicit with the aid of (1B.4j) . Eq. (IB. 71) becomes 
condition (12. 3p of the main text. 
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